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Scientists have been considering the Kuramoto model to understand the mechanism behind the 
appearance of collective behaviour, such as frequency synchronisation (FS) as a paradigm, in real- 
world networks with a finite number of oscillators. A major current challenge is to obtain an 
analytical solution for the phase-angles. Here, we provide an approximate analytical solution for 
this problem by deriving a master solution for the finite-size Kuramoto model, without imposing 
any restriction on the distribution of the natural frequencies of the oscillators. The master solution 
embodies all particular solutions of the finite-size Kuramoto model for any frequency distribution 
and coupling strength larger than the critical one. Furthermore, we present a criterion to determine 
the stability of the FS solution. This allows one to analytically infer the relationship between the 
physical parameters and the stable behaviour of networks. 

PACS numbers: 89.75.-k, 05.45.Xt, 05.45.-a 


Networks of coupled oscillators provide a pragmatic 
model to describe basic behaviour of natural and tech¬ 
nological systems, such as biological networks , so¬ 
cial networks Q, computer networks ii and power 
grids @,0- A significant phenomenon emerging from 
coupled oscillators is the synchronisation of the oscilla¬ 
tors’ rhythms [Mil- Since 1958, when Norbert Wiener 
[ 1 ^ suggested the presence of collective behaviour of neu¬ 
rons in the brain, finding the mechanism and the con¬ 
ditions for the appearance of synchronisation in com¬ 
plex networks has attracted attention of many scientists. 
Kuramoto (l3l - [T5l| proposed a mathematically tractable 
model by considering a network of phase oscillators cou¬ 
pled with an all-to-all topology through a sine function, 
with each oscillator possessing a constant natural fre¬ 
quency. Recently, scientists devoted great effort to un¬ 
derstand the collective behaviour and synchronisation of 
the Kuramoto model and its connections to natural sys¬ 
tems [l6l - [^ . 

However, most of the existing analytical results have 
considered networks with an infinite number of oscillators 
and have imposed restrictions on the distributions of nat¬ 
ural frequencies of oscillators. To the best of our knowl¬ 
edge, no analytical work has been proposed that allows 
one to determine the synchronous solution of the phase- 
angles for the finite-size Kuramoto model without restric¬ 
tion on natural frequencies. Although the conditions 
for the set of natural frequencies that provides stable 
synchronisation is obtained in [ 2 ^, the solution for the 
phase-angles is not given. But, solving the phase-angles 
is of fundamental significance in real-world systems. For 
example, in power grids [a, which can be described 
by a Kuramoto-like model [24| , one is not only interested 
in understanding under which conditions the frequency 
synchronisation (FS) among oscillators emerges, but also 
in knowing the phase-angles after synchronisation. The 
phase-angles are important variables for monitoring gen¬ 
erators and developing control strategies for the power 


grid . Similarly, in the research of Josephson junctions 
[^ . where a Kuramoto-Sakaguchi model is consid¬ 
ered, the phase-angles of a synchronous state contain the 
information of the wave-function phase difference across 
every Josephson junction. 

In this work, we present an analytical method to ap¬ 
proximately calculate the phase-angles for the finite-size 
Kuramoto model when the frequencies are synchronised. 
Compared to other works, our method does not require 
any restriction on the distribution of the natural frequen¬ 
cies of oscillators. Our solution, shown in Eq. dH, links 
directly the FS solution and the physical parameters of 
the network. Remarkably, the solution is independent 
of the network size, and only depends on the natural 
frequencies and the coupling strength. In addition, we 
provide an approximate criterion, shown in Eq. (1261) . to 
analytically predict whether a finite number of oscilla¬ 
tors are able to emerge into a stable FS, even without 
knowing the FS solution explicitly. 

We use In (On) to denote the TV x 1 vector with all 
elements equal to one (zero), In to indicate the index set 
{1,2,-- - ,TV}, and M to represent the set of real numbers. 
Given a vector a with TV elements, we use a = ^ 
to denote the mean value of a. 

The finite-size Kuramoto, describing the dynamical be¬ 
haviour of TV phase oscillators in an all-to-all network, is 
given by the equation, 

K ^ 

Oi = n^ +—'^sin(ej -Qi), ieiN, ( 1 ) 

T=i 

where TV is a finite positive integer number, K is 
the coupling strength, and the TV x 1 vectors, fl = 
[Hi, H 2 , ••• , ^Nf and 0 = [ 01 , 02 , ••• , Sn]^, 
denote the natural frequencies and instantaneous phases 
of the oscillators, respectively. We rewrite Eq. m in a 
rotating frame by letting 6 = Q — and d; = Q — 
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such that a; = 0. Thus, Eq. o becomes 


9^ = uJi +—'^sm{6j - 6i), i&lN- ( 2 ) 


N 


N 




The order parameter r € [0,1] is defined by [3 [3 
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i=i 


where r and V' are calculated by equating the real and 
imaginary parts in Eq. ®, namely, 
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^sin6lj 


and 


i=i 


tp = arctan 


y;^iCos6ij- 


\ ^ 

1 \ 
-^cos^J , (4) 


(5) 


We note that ip is not the average phase in a network 
with finite oscillators (i.e., ip ^ 6). This is different from 
the situation with infinite number of oscillators HI] , 
where ip = 9 is the global mean field. 

Multiplying to both sides of Eq. m and equating 
the imaginary part, we have 


r sin {ip — 9i) 


1 ^ 

— ^ sin (0^-- 6»,), 


( 6 ) 


which results in 


9^ = uji +Krs\n{ip - 9i), iGIn- (7) 
The oscillators described by Eq. ([21) emerge into ES if 

M- 


0i — = 0 as t —>• oo, Vi, j € In- (8) 

We have, from Eq. (I2|), that 

Since w = 0, to satisfy Eq. ([8|) we require that 9 = InoJ = 
On- Thus, in the rotating framework, the finite-size Ku- 
ramoto model in ES can be described by 


N 


4 = Wi -I- — ^ sin(0j - 9i), 
^ 1=1 

9i =0, Vi G In- 


(9) 


where u is calculated from y 1 — (ujj/u)'^ = 

The stability of the ES solution of the finite-size Ku- 
ramoto model can be studied by the Lyapunov function, 

E = (11) 

where the time derivative of E along the trajectories of 
Eq. dlj) is E = -f Ef^i<js;ArCOs(0* - 9j){9i - 9jf. A 
sufficient condition for the stability of the ES solution is 
that E < 0, implying cos(9i — 9j) > 0, Vi, j G In- This 
means that \9i — 9j \ < Vi, j G In- We denote the sta¬ 
ble critical coupling strength of the finite-size Kuramoto 
model by Ks defined as 

Ks ■-= min{iG : K ^ Kc, \9, - 9,\ < J, Vi, j G In}- 

( 12 ) 

Thus, as i —>■ oo, the oscillators are attracted to the stable 
ES if AT ^ Ks- We define Ks ^ Kq, since K = Kq only 
ensures the existence of the ES solution Q , but K ^ Ks 
provides a condition for its stability based on the Lya- 
punov function. The existence of Ks has been studied in 
m, 0 and [13. 

Given K = K* ^ Ks, assume 0® is a stable so¬ 
lution of the finite-size Kuramoto model in FS. Then, 
0« = 0® + 1^^, V^ G M is also a stable solution cor¬ 
responding to AT = AT*, since the phase differences are 
independent of ^ . This means that there are an infinite 
number of stable solutions for the finite-size Kuramoto 
model in a certain FS state. We use , correspond¬ 
ing to AT = AT*, to denote the infinite-dimensional stable 
solution space. Actually, 0 = In^, V^ G K, is the ho¬ 
mogeneous solution of Eq. obtained by setting the 
non-homogeneous terms to be zero, i.e., uj = On- 
and 0^ are particular solutions of the non-homogeneous 
Eq. (I2|). Our goal is to find one of the particular solutions, 
which can be analytically expressed, such that Eq. m is 
analytically solvable. 

We define the master solution of Eq. ([9]) as 

^*=0«-lw^«GSf., (13) 

where ip^ is calculated from Eq. m as 0 = 0^. There are 
three characteristics for which we call cp* as the master 
solution: (i) cp* is analytically expressible; (ii) (p* is iden¬ 
tical V^ G K, i.e., (p* is independent of (iii) ip* = 0, 
where ip* is calculated from Eq. as 9 = (p*. Next, we 
will prove the three characteristics of cp*. 

Considering Eq. ([6]), Eq. ([9]) can be transformed into 


In general, the finite-size Kuramoto model reaches FS if 
the coupling strength is larrar than a critical value, A'q 
M- Verwoerd and Mason 0 provided an algorithm to 
exactly calculate Kc from 


Kc = 


1 

N 



( 10 ) 


= K*r* sia(P*, Vi G In, (14) 

where r* is calculated by multiplying on both sides 

of Eq. ([3]), namely, 

AT JV 

r* = cos(<(>*). (15) 

i=i i=i 
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We define > where [ft^tnax) is 

the minimum (maximum) (jf. We have two observations 
about T^,: firstly, its length |T^. | < § due to 4>* S ; 
secondly, 0 € , because Za = 0 implies uJmin ^ 0 

and Umax > 0, resulting in ^ 0 and > 0 

from Eq. (HI- These two characteristics of T^, restrict 
T^. C [—f, f]- Thus, we get the analytical expression 
of the master solution ((p*) from Eq. (iMl) . 

ip* = arcsin-^^, Vi € In, (16) 

which is independent of thus, (i) and (ii) follows. 
From Eq. ((H]) we have sin<(>)‘ = = 6- 

Substituting this into Eq. m, we have tp* = 0, which 
proves (iii). Our goal is to find the master solution 
{(p*) as K = K* ^ Ks for the finite-size Kuramoto 
model. We know T^. C [—f, f], implying cos cpj ^ 
0, Vj G. In- Thus, the order parameter in Eq. (US is 
calculated as r* = ^ \J^ ~ sin(((>*)^. Considering 

sin^* = Uj/{K*r*) obtained from Eq. ([14]), r* can be 
expressed by a transcendental equation, namely. 



In order to obtain an approximate solution for r*, we 
construct a new model for the original system. We rela¬ 
bel u such that ui Vi U 2 ^ ^ un, and split u into 

two groups. One group is u' = [ui,U 2 , ■ ■ ■ jUn']'^, where 
N' = ^ {N' = ■^^^) if N is even (odd), and the other 
group is u” = [un'+i,'^n'+ 2 ,' ■ ■ ,'^n]’^- We note that, 

j=l ^3 + 2Zj=N' + l ^3 = Si=l ^3= ^ 

overlineu" Vi 0 when N is even, and u' ~ —u" ^ 0 
when N is odd. For simplicity, we indistinctly denote 
w' ~ —u" for both cases. 

When all of the oscillators emerge into stable FS, our 
model treats the whole system as two oscillators in stable 
FS. The natural frequencies of the two oscillators are u' 
and u". The two oscillators also follow the original FS 
Kuramoto model equations, namely, 


u' = K*r*' sin^'. 

(18) 

^ = /CV*'sin^", 

(19) 


where r*' is the order parameter for the two oscilla¬ 
tors. The two oscillators are in stable FS, thus we have 
(p', cp" € [~f j f]: which is obtained from the analysis of 
the Lyapunov function [Eq. (ITT]) ] for the two oscillators. 
Then the order parameter is r*' = 1/2 (cos^'-i-cos0") = 

1/2^1 - p'/{K*r*)Y -I- 1/2^1 - p"/{K*r*)Y . Further 
considering |a;^| ~ \uj"\ we have 



whose solution is 




Ai = 


V2, 


A 


i + \ 1- 


4a;' 


K*2 


K* > 2\u'\ 


( 21 ) 




\ 


1 - \ 1 - 


4a;' 


K*2 


K* ^ 2\u'\, (22) 


where we have A 1 A 2 = ~ Substituting this 

condition into Eq. (m and Eq. (HI, we get — A 1 A 2 = 
r*' sincp' and A 1 A 2 = r*'sin^". If r*' ss A 2 , we have (p' ~ 
— arcsin(Ai) and (p" ~ arcsin(Ai). Because ^ ^ Ai < 1, 
we approximately have — f ^ ^ — f and ^ 

This means \(p' — (p"\ ^ and Ai grows larger as K in¬ 
creases from Ks resulting in a growth oi\(p' — (p"\. How¬ 


ever, \(p' — (p"\ 




implies instability of the FS solution 


of the two oscillators, which can be understood by the 
Lyapunov function in Eq. (HD for the two oscillators. 
This means that r*' Ri A 2 describes an unstable FS solu¬ 
tion. On the other hand, r*' k, Ai ensures the stability 
of the FS solution. Thus, we let r* Ri r*' Ki Ai be the 
approximation for the order parameter in Eq. dnD, then 
the approximation {(p**) for the master solution {(p*) in 
Eq. (|T^ is 


(p** = arcsin Vi £ In (23) 

Consequently, the approximate stable solution of the 
finite-size Kuramoto model in FS is 

9i ss arcsin —h Vi G In, G R. (24) 

i\ Ai 

We use e with element Ci = \(p** — (p*\, Vi G In, to de¬ 
note the absolute error between the approximate master 
solution [(p** in Eq. d^ ] and the numerical one [cp* in 
Eq. (ITBI) ]. and a to denote the standard deviation of e, 

defined as cr = To demonstrate the 

effectiveness of our method for a network with four os¬ 
cillators, we show numerical results in Fig. |T] is large 
Vi G In, when K < Ks [(p** ^ (p* in Fig. [1] (a)], and 
small when K = Ks [$** ~ (p* in Fig. |T] (b)]. The ap¬ 
proximate order parameter, Ai, is close to the numerical 
one as K ^ Ks [Fig. |T] (c)]. The average absolute error 
and standard deviation of the absolute error are small 
(e < 10“^, cr < 10“^) when K ^ Ks, and decreases 
rapidly as K is increased [Fig. [1] (d)]. This means that 
our method can calculate the master solution almost ex¬ 
actly in an analytical way for this network. 

Figured] indicates that our method works well for net¬ 
works with different Q distributions (random, normal and 
exponential) and different number of oscillators {N in¬ 
creasing from 3 to 200). We observe that e and cr are 
larger when N is smaller, but they decrease quickly as N 
increases. Furthermore, e and cr in Fig. d) are obtained 
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FIG. 1. [Colour online] Results for a network of 4 oscil¬ 
lators with the vector of natural frequency given by = 
[—2, —1,1, 2]^. (a) K < Ks, the red triangles and blue circles 
indicate the approximate master solution ((/>**) in Eq. (1231) 
and the numerical one (0*) in Eq. (11611 . respectively, (b) 
K = Ks, </>** indicated by red triangles is close to 0* indi¬ 
cated by blue circles, (c) The change of the order parame¬ 
ter and its approximation with respect to K. The red solid 
line indicates the numerical result of the order parameter in 
Eq. (1171) . (an average value of results from 2000 simulations 
with different initial phases). The blue solid line and green 
solid line indicate the change of Ai in Eq. ini) and A 2 in 
Eq. (1221) . respectively, as K ^ 2|a;'|, where 2 |cj'| is the lowest 
bound of K in Eq. 1211) and Eq. (1221) . The magenta (cyan, 
black) dash line represents 2\uj'\ {Kc, Ks)- (d) The change of 
the average of the absolute error (red solid line) and standard 
deviation of the absolute error (blue solid line) between (j)** 
and (j)* as a function of K when K ^ 2|a;'|. 


at K = Ks, meaning that they are the largest value ob¬ 
tained for each simulation. In other words, smaller e and 
cr can be obtained if we increase K for any given N and 
any Q distribution, since larger K implies smaller abso¬ 
lute error between the approximate order parameter and 
the numerical one as shown in Fig. [T] (c), which further 
implies smaller absolute error between (j)** and (jf. 

A sufficient condition to ensure the stability of the FS 
solution was proposed by @ as 

K > = (25) 

r 

Taking Ai in Eq. (ED as the approximation of r, we get 
a sufficient condition for the stability of the FS solution, 
namely. 


K> Ka = Vf G In- (26) 

Ai 

Equation (1261) is useful to approximately forecast whether 
the system is able to get into a stable FS state in an ana¬ 
lytical way, without solving the differential equations. In 



FIG. 2. [Golour online] Study of the average absolute error e 
in (a), and the standard deviation a in (b) between the ap¬ 
proximate master solution (0**) and the numerical one {(!>*) 
at K = Ks- N increases from 3 to 200. The red line (green 
line, black line) corresponds to the normal distribution (ran¬ 
dom distribution, exponential distribution) of H. Results are 
based on the average value of results from 100 simulations for 
each distribution. 


other words, it reveals the relationship between the phys¬ 
ical parameters (average natural frequency and coupling 
strength) and the stable behaviour of oscillator networks 
described by the finite-size Kuramoto model explicitly. 





FIG. 3. [Golour online] The stable critical coupling strength 
of the finite-size Kuramoto model. N increases from 3 to 200. 
Ks is calculated by Eq. (O, Kp hy Eq. (1^ . Ka by Eq. TIEh - 
(a) Q follows a normal distribution, (b) Q follows a random 
distribution, (c) Q follows an exponential distribution. The 
figures are drawn based on average results of 30 simulations 
for each distribution. 

Figure [3] shows the effectiveness of our condition to 
gauge the stability for the FS solution of the finite-size 
Kuramoto model, considering different 12 distributions 
and different number of oscillators. If 12 obeys a normal 
or exponential distribution [Figs. [3] (a) and[3](c), respec¬ 
tively] , Ka in Eq. (1261) coincides remarkably with Kp in 
Eq. (1251) . and Ka > Ks- If 12 obeys a random distribu¬ 
tion [Fig.[3](b)], Ka in Eq. (1251) is close to Ks in Eq. (1T21) . 
This means Eq. (1261) is effective to approximately deter¬ 
mine the stability for the FS solution of the finite-size 
Kuramoto model. 

In this paper, we studied the finite-size Kuramoto 
model [Eq. (I2D], including its analytical solution for fre¬ 
quency synchronisation (FS) and its stable behaviour. 
Our approximate FS solution takes a simple form 
[Eq. (1231) ]. which, surprisingly, is independent of the net¬ 
work size. These significant results is a consequence of a 
mathematical insight expressed in Eq. dHD and a physi¬ 
cal insight in the model leading to Eq. (1T21) and Eq. (1191) . 
Among an infinite number of FS solutions, we have un- 





















5 


derstood that there is a particular one, the master solu¬ 
tion in Eq. which allows one to calculate all of the 
others. We have also understood that the FS in the finite- 
size Kuramoto model is approximately characterised by 2 
clusters of oscillators that can be effectively described as 


2 coupled oscillators in Eq. (1181) and Eq. dm. Further¬ 
more, we developed a condition to approximately predict 
the stability for the FS solution of the finite-size Ku¬ 
ramoto model in an analytical way. This condition allows 
one to easily infer the relationship between the physical 
parameters and the stable behaviour of networks. 
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